#ifndef PCPS_ACCUMULATOR_PROPERTIES_HEADER
#define PCPS_ACCUMULATOR_PROPERTIES_HEADER

#include <src/pcps.h>
#include <src/accumulator.h>
#include <src/operator.h>
#include <src/local_correlation.h>
#include <src/icosidodecahedron.h>

namespace pcps {

  // forward declaration of thread accumulator
  template <class SCALAR, class PART, class REF> class PropertiesThreadAccum;

  //--------------------------------------------------------------------------------
  // pcps::PropertiesProcessAccum -- Accumulates correlation functions.
  //--------------------------------------------------------------------------------

  template <class SCALAR, class PART, class REF> class PropertiesProcessAccum : public pcps::ProcessAccum<SCALAR, PART, REF> {

    private:

      // data members
      const int _nsites;                        // number of sites in the lattice
      const int _configs_per_site;              // the number of configurations in a single site's Hilbert space
      const bool _spin_pair_correlations;       // whether to compute pair spin correlations for all pairs
      const bool _ssss_correlations;            // whether to compute the four site Si.Sj Sk.Sl correlations
      const bool _spin_triple_products;         // whether to compute icosidodecahedron spin triple product correlations
      const bool _spin_triangle_3sz;            // whether to compute icosidodecahedron sz*sz*sz correlations for each triangle
      std::vector<SCALAR> _process_sz;          // array to hold this process's contribution to < Sz_i > for each site i
      std::vector<SCALAR> _process_sz_dot_sz;   // array to hold this process's contribution to < Sz_i Sz_j > for each pair of sites i, j
      std::vector<SCALAR> _process_s_dot_s;     // array to hold this process's contribution to < S_i . S_j > for each pair of sites i, j
      std::vector<SCALAR> _process_ssss;        // array to hold this process's contribution to < S_i.S_j S_k.S_l >
      std::vector<SCALAR> _process_stp;         // array to hold this process's contribution to < S_i.(S_j x S_k) > for each icosidodecahedron triangle
      std::vector<SCALAR> _process_3sz;         // array to hold this process's contribution to < Sz_i * Sz_j * Sz_k > for each icosidodecahedron triangle
      std::vector<std::vector<int> > _triangles;// vectors holding the triangles' sites for the icosidodecahedron
      std::set<int> _ico_a_sites;               // set of a sublattice sites on the icosidodecahedron
      std::set<int> _ico_b_sites;               // set of b sublattice sites on the icosidodecahedron
      std::set<int> _ico_c_sites;               // set of c sublattice sites on the icosidodecahedron
      const pcps::action_value_object<SCALAR, pcps::OpRL<SCALAR>, PART> _avoRL;  // object for computing numerical coefficient
                                                                                  // for action of raising/lowering operator

    public:

      // give the corresponding thread accumulator access to private and protected members
      friend class pcps::PropertiesThreadAccum<SCALAR, PART, REF>;

      // constructor
      PropertiesProcessAccum(const pcps::Input & userinp, const pcps::Wavefunction<SCALAR> & wfn, const long int sample_length)
        : pcps::ProcessAccum<SCALAR, PART, REF>(userinp, sample_length),
          _nsites(userinp.nsites()),
          _configs_per_site(userinp.configs_per_site()),
          _spin_pair_correlations( userinp.spin_pair_correlations() ),
          _ssss_correlations( userinp.ssss_correlations() ),
          _spin_triple_products( userinp.spin_triple_products() ),
          _spin_triangle_3sz( userinp.spin_triangle_3sz() ),
          _avoRL(userinp.configs_per_site())
      {

        if ( _spin_triple_products && userinp.lattice_type() != pcps::keyword::Lattice_Icosidodecahedron )
          throw pcps::Exception("spin triple products are only defined on the icosidodecahedron lattice");

        if ( _spin_triangle_3sz && userinp.lattice_type() != pcps::keyword::Lattice_Icosidodecahedron )
          throw pcps::Exception("spin triangle sz*sz*sz products are only defined on the icosidodecahedron lattice");

        if ( _ssss_correlations && !_spin_pair_correlations )
          throw pcps::Exception("\"ssss_correlations = true\" requires that \"spin_pair_correlations = true\"");

        // initialize arrays
        const SCALAR z = pcps::zero<SCALAR>();
        if ( _spin_pair_correlations ) {
          _process_sz.assign(         size_t(_nsites)                         , z );
          _process_sz_dot_sz.assign(  size_t(_nsites*_nsites)                 , z );
          _process_s_dot_s.assign(    size_t(_nsites*_nsites)                 , z );
        }
        if ( _ssss_correlations ) {
          _process_ssss.assign(       size_t(_nsites*_nsites*_nsites*_nsites) , z );
        }
        if ( _spin_triple_products ) {
          _process_stp.assign(        size_t(20)                              , z );
        }
        if ( _spin_triangle_3sz ) {
          _process_3sz.assign(     size_t(20)                              , z );
        }

        // get triangles
        if ( _spin_triple_products || _spin_triangle_3sz ) {
          std::set<std::set<int> > tri_sets = pcps::icosidod_get_three_pent_triples();
          _triangles.resize(tri_sets.size());
          std::set<std::set<int> >::const_iterator s = tri_sets.begin();
          std::vector<std::vector<int> >::iterator v = _triangles.begin();
          for ( ; v != _triangles.end() && s != tri_sets.end(); v++, s++)
            v->assign(s->begin(), s->end());
          _ico_a_sites = pcps::icosidod_get_a_sites();
          _ico_b_sites = pcps::icosidod_get_b_sites();
          _ico_c_sites = pcps::icosidod_get_c_sites();
        }

        this->reset();

      }

      // function to get an icosidodecahedron triangle's sites in a,b,c order
      void abc_order_ico_triangle_sites(const int s0, const int s1, const int s2, int & a, int & b, int & c) const {

        a = -1;
        b = -1;
        c = -1;

        if      ( _ico_a_sites.count(s0) > 0 ) a = s0;
        else if ( _ico_b_sites.count(s0) > 0 ) b = s0;
        else if ( _ico_c_sites.count(s0) > 0 ) c = s0;

        if      ( _ico_a_sites.count(s1) > 0 ) a = s1;
        else if ( _ico_b_sites.count(s1) > 0 ) b = s1;
        else if ( _ico_c_sites.count(s1) > 0 ) c = s1;

        if      ( _ico_a_sites.count(s2) > 0 ) a = s2;
        else if ( _ico_b_sites.count(s2) > 0 ) b = s2;
        else if ( _ico_c_sites.count(s2) > 0 ) c = s2;

        if ( a < 0) throw pcps::Exception("failed to find the A sublattice site in abc_order_ico_triangle_sites");
        if ( b < 0) throw pcps::Exception("failed to find the B sublattice site in abc_order_ico_triangle_sites");
        if ( c < 0) throw pcps::Exception("failed to find the C sublattice site in abc_order_ico_triangle_sites");

      }

      // function to create a thread accumulator linked to this process accumulator
      pcps::ThreadAccum<SCALAR, PART, REF> * create_thread_accum(const int tid, pcps::Arena & arena);

      // function to reset the accumulator
      void reset() {
        if ( _spin_pair_correlations ) {
          pcps::xscal(                        _nsites, pcps::zero<SCALAR>(), &_process_sz.at(0),        1);
          pcps::xscal(                _nsites*_nsites, pcps::zero<SCALAR>(), &_process_sz_dot_sz.at(0), 1);
          pcps::xscal(                _nsites*_nsites, pcps::zero<SCALAR>(), &_process_s_dot_s.at(0),   1);
        }
        if ( _ssss_correlations ) {
          pcps::xscal(_nsites*_nsites*_nsites*_nsites, pcps::zero<SCALAR>(), &_process_ssss.at(0),      1);
        }
        if ( _spin_triple_products ) {
          pcps::xscal(                             20, pcps::zero<SCALAR>(), &_process_stp.at(0),       1);
        }
        if ( _spin_triangle_3sz ) {
          pcps::xscal(                             20, pcps::zero<SCALAR>(), &_process_3sz.at(0),       1);
        }
      }

      // function to print information during stochastic iterations
      void print_iter_info(const pcps::Input & userinp,
                           pcps::Wavefunction<SCALAR> & wfn,
                           const int iter,
                           std::vector<boost::shared_ptr<pcps::ProcessAccum<SCALAR, PART, REF> > > & accumulators) {

        if (this->myrank() == 0) {
          const std::string filename = this->get_file_name(iter);
          std::cout << boost::format("properties accumulator:     properties written to file \"%s\"") % filename << std::endl;
          this->write_to_file(filename);
        }

      }

      // function to get the filename to write the properties to
      std::string get_file_name(const int iter) const {
        std::string retval = (boost::format("%i") % iter).str();
        if (retval.size() < 8)
          retval.insert(retval.begin(), 8 - retval.size(), '0');
        retval = (boost::format("properties_iter_%s.txt") % retval).str();
        return retval;
      }

      // function to write properties to a file
      void write_to_file(const std::string & filename) const {

        // open the output file
        std::ofstream outfile(filename.c_str());
        if (!outfile.is_open())
          throw pcps::Exception( (boost::format("failed to open file \"%s\" for writing properties") % filename).str() );

        if ( _spin_pair_correlations ) {

          // write the Sz averages
          outfile << boost::format("%6s  %40s") % "site i" % "< Sz_i >" << std::endl;
          for (int i = 0; i < _nsites; i++) {
            outfile << boost::format("%6i  %40s")
                       % i
                       % pcps::format_number("%14.6f", _process_sz[i] )
                    << std::endl;
          }
          outfile << std::endl;

          // write the Sz_i Sz_j averages
          outfile << boost::format("%6s  %6s  %40s") % "site i" % "site j" % "< Sz_i Sz_j >" << std::endl;
          for (int i = 0; i < _nsites; i++)
          for (int j = 0; j < _nsites; j++) {
            outfile << boost::format("%6i  %6i  %40s")
                       % i
                       % j
                       % pcps::format_number("%14.6f", _process_sz_dot_sz[i*_nsites+j] )
                    << std::endl;
          }
          outfile << std::endl;

          // write the S_i.S_j averages
          outfile << boost::format("%6s  %6s  %40s") % "site i" % "site j" % "< S_i . S_j >" << std::endl;
          for (int i = 0; i < _nsites; i++)
          for (int j = 0; j < _nsites; j++) {
            outfile << boost::format("%6i  %6i  %40s")
                       % i
                       % j
                       % pcps::format_number("%14.6f", _process_s_dot_s[i*_nsites+j] )
                    << std::endl;
          }
          outfile << std::endl;

        }

        if ( _ssss_correlations ) {

          // write the S_i.S_j S_k.S_l averages
          outfile << boost::format("%6s  %6s  %6s  %6s  %40s") % "site i" % "site j" % "site k" % "site l" % "< ( S_i . S_j ) ( S_k . S_l ) >" << std::endl;
          for (int i = 0; i < _nsites; i++)
          for (int j = 0; j < _nsites; j++)
          for (int k = 0; k < _nsites; k++)
          for (int l = 0; l < _nsites; l++) {
            outfile << boost::format("%6i  %6i  %6i  %6i  %40s")
                       % i
                       % j
                       % k
                       % l
                       % pcps::format_number("%14.6f", _process_ssss[_nsites*(_nsites*(_nsites*i+j)+k)+l] )
                    << std::endl;
          }
          outfile << std::endl;

        }

        if ( _spin_triple_products ) {

          outfile << boost::format("spin triangle triple products:") << std::endl;
          for (int t = 0; t < _triangles.size(); t++) {
            int i, j, k;
            this->abc_order_ico_triangle_sites(_triangles.at(t).at(0), _triangles.at(t).at(1), _triangles.at(t).at(2), i, j, k);
            outfile << boost::format(" %6i %6i %6i %55s") % i % j % k % pcps::format_number("%20.12f", _process_stp.at(t)) << std::endl;
          }
          outfile << std::endl;

        }

        if ( _spin_triangle_3sz ) {

          outfile << boost::format("spin triangle sz*sz*sz:") << std::endl;
          for (int t = 0; t < _triangles.size(); t++) {
            const int i = _triangles.at(t).at(0);
            const int j = _triangles.at(t).at(1);
            const int k = _triangles.at(t).at(2);
            outfile << boost::format(" %6i %6i %6i %55s") % i % j % k % pcps::format_number("%20.12f", _process_3sz.at(t)) << std::endl;
          }
          outfile << std::endl;

        }

        // close the file
        outfile.close();

      }

  };

  //--------------------------------------------------------------------------------
  // pcps::PropertiesThreadAccum -- Accumulates correlation functions.
  //--------------------------------------------------------------------------------

  template <class SCALAR, class PART, class REF> class PropertiesThreadAccum : public pcps::ThreadAccum<SCALAR, PART, REF> {

    private:

      // pointer to the corresponding process accumulator
      pcps::PropertiesProcessAccum<SCALAR, PART, REF> * const _ppacc;

      // data members
      const int _nsites;            // number of sites in the lattice
      const int _configs_per_site;  // the number of configurations in a single site's Hilbert space
      SCALAR * _thread_sz;          // array to hold this thread's contribution to < Sz_i > for each site i
      SCALAR * _thread_sz_dot_sz;   // array to hold this thread's contribution to < Sz_i Sz_j > for each pair of sites i, j
      SCALAR * _thread_s_dot_s;     // array to hold this thread's contribution to < S_i.S_j > for each pair of sites i, j
      SCALAR * _thread_ssss;        // array to hold this thread's contribution to < S_i.S_j S_k.S_l>
      SCALAR * _thread_stp;         // array to hold this thread's contribution to < S_i.(S_j x S_k) > for each icosidodecahedron triangle
      SCALAR * _thread_3sz;         // array to hold this thread's contribution to < Sz_i * Sz_j * Sz_k > for each icosidodecahedron triangle
      int * _workspace;             // array for use when computing local correlation contributions

    public:

      // constructor
      PropertiesThreadAccum(pcps::PropertiesProcessAccum<SCALAR, PART, REF> * const ppacc, const int tid, pcps::Arena & arena)
        : pcps::ThreadAccum<SCALAR, PART, REF>(ppacc, tid),
          _ppacc(ppacc),
          _nsites(ppacc->_nsites),
          _configs_per_site(ppacc->_configs_per_site)
      {

        // initialize arrays
        const SCALAR z = pcps::zero<SCALAR>();
        if ( this->_ppacc->_spin_pair_correlations ) {
          _thread_sz         = arena.allocate_array<  SCALAR>( size_t(_nsites)                         , z );
          _thread_sz_dot_sz  = arena.allocate_array<  SCALAR>( size_t(_nsites*_nsites)                 , z );
          _thread_s_dot_s    = arena.allocate_array<  SCALAR>( size_t(_nsites*_nsites)                 , z );
        }
        if ( this->_ppacc->_ssss_correlations ) {
          _thread_ssss       = arena.allocate_array<  SCALAR>( size_t(_nsites*_nsites*_nsites*_nsites) , z );
        }
        if ( this->_ppacc->_spin_triple_products ) {
          _thread_stp        = arena.allocate_array<  SCALAR>( size_t(20)                              , z );
        }
        if ( this->_ppacc->_spin_triangle_3sz ) {
          _thread_3sz        = arena.allocate_array<  SCALAR>( size_t(20)                              , z );
        }
        _workspace         = arena.allocate_array<     int>( size_t(_nsites)                         , 0 );

        // reset data
        this->reset();

      }

      // function to reset the accumulator
      void reset() {
        if ( this->_ppacc->_spin_pair_correlations ) {
          pcps::xscal(                        _nsites, pcps::zero<SCALAR>(), _thread_sz,        1);
          pcps::xscal(                _nsites*_nsites, pcps::zero<SCALAR>(), _thread_sz_dot_sz, 1);
          pcps::xscal(                _nsites*_nsites, pcps::zero<SCALAR>(), _thread_s_dot_s,   1);
        }
        if ( this->_ppacc->_ssss_correlations ) {
          pcps::xscal(_nsites*_nsites*_nsites*_nsites, pcps::zero<SCALAR>(), _thread_ssss,      1);
        }
        if ( this->_ppacc->_spin_triple_products ) {
          pcps::xscal(                             20, pcps::zero<SCALAR>(), _thread_stp,       1);
        }
        if ( this->_ppacc->_spin_triangle_3sz ) {
          pcps::xscal(                             20, pcps::zero<SCALAR>(), _thread_3sz,       1);
        }
      }

      // function to accumulate the correlation functions of the old configuration
      void accumulate(const long int iter, 
                      const long int count,
                      const bool accept,
                      const bool finished,
                      const SCALAR coeff_ratio,
                      const int * const new_config,
                      const int * const old_config,
                      pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf,
                      pcps::ThreadAccum<SCALAR, PART, REF> ** accumulators,
                      const int n_accumulators) {

        if (accept || finished) {

          // compute correlation functions
          if ( this->_ppacc->_spin_pair_correlations ) {
            SCALAR * const ssss = ( this->_ppacc->_ssss_correlations ? _thread_ssss : NULL );
            pcps::local_correlation_functions<SCALAR, PART, REF>(_nsites, _configs_per_site, count, old_config, _workspace, _ppacc->_avoRL, cvf,
                                                                 _thread_sz, _thread_sz_dot_sz, _thread_s_dot_s, ssss);
          }

          // compute spin triple product correlations for the icosidodecahedron's triangles
          if ( this->_ppacc->_spin_triple_products ) {
            this->accumulate_spin_triple_products(count, old_config, cvf);
          }

          // compute sz*sz*sz for the icosidodecahedron's triangles
          if ( this->_ppacc->_spin_triangle_3sz ) {
            this->accumulate_spin_triangle_3sz(count, old_config, cvf);
          }

        }

        if (finished) {

          // add this thread's contribution to the process's total
          pthread_mutex_lock(this->_pacc->mutex());
          const SCALAR total_samples = pcps::unity<SCALAR>() * this->_pacc->total_samples();
          if ( this->_ppacc->_spin_pair_correlations ) {
            pcps::xaxpy(                        _nsites, 1.0 / total_samples, _thread_sz,        1, &_ppacc->_process_sz.at(0),        1);
            pcps::xaxpy(                _nsites*_nsites, 1.0 / total_samples, _thread_sz_dot_sz, 1, &_ppacc->_process_sz_dot_sz.at(0), 1);
            pcps::xaxpy(                _nsites*_nsites, 1.0 / total_samples, _thread_s_dot_s,   1, &_ppacc->_process_s_dot_s.at(0),   1);
          }
          if ( this->_ppacc->_ssss_correlations ) {
            pcps::xaxpy(_nsites*_nsites*_nsites*_nsites, 1.0 / total_samples, _thread_ssss,      1, &_ppacc->_process_ssss.at(0),      1);
          }
          if ( this->_ppacc->_spin_triple_products ) {
            pcps::xaxpy(                             20, 1.0 / total_samples, _thread_stp,       1, &_ppacc->_process_stp.at(0),       1);
          }
          if ( this->_ppacc->_spin_triangle_3sz ) {
            pcps::xaxpy(                             20, 1.0 / total_samples, _thread_3sz,       1, &_ppacc->_process_3sz.at(0),       1);
          }
          pthread_mutex_unlock(this->_pacc->mutex());

          // wait for other threads
          pcps::pthread_barrier_wait(this->_pacc->barrier());

          // combine results from all processors on root process
          if (this->_tid == 0) {

            const MPI::Comm & comm = MPI::COMM_WORLD;

            if ( this->_ppacc->_spin_pair_correlations ) {
              pcps::xcopy(                        _nsites, &_ppacc->_process_sz.at(0),        1, _thread_sz,        1);
              pcps::xcopy(                _nsites*_nsites, &_ppacc->_process_sz_dot_sz.at(0), 1, _thread_sz_dot_sz, 1);
              pcps::xcopy(                _nsites*_nsites, &_ppacc->_process_s_dot_s.at(0),   1, _thread_s_dot_s,   1);
            }
            if ( this->_ppacc->_ssss_correlations ) {
              pcps::xcopy(_nsites*_nsites*_nsites*_nsites, &_ppacc->_process_ssss.at(0),      1, _thread_ssss,      1);
            }
            if ( this->_ppacc->_spin_triple_products ) {
              pcps::xcopy(                             20, &_ppacc->_process_stp.at(0),       1, _thread_stp,       1);
            }
            if ( this->_ppacc->_spin_triangle_3sz ) {
              pcps::xcopy(                             20, &_ppacc->_process_3sz.at(0),       1, _thread_3sz,       1);
            }

            if ( this->_ppacc->_spin_pair_correlations ) {
              pcps::xscal(                        _nsites, pcps::zero<SCALAR>(), &_ppacc->_process_sz.at(0),        1);
              pcps::xscal(                _nsites*_nsites, pcps::zero<SCALAR>(), &_ppacc->_process_sz_dot_sz.at(0), 1);
              pcps::xscal(                _nsites*_nsites, pcps::zero<SCALAR>(), &_ppacc->_process_s_dot_s.at(0),   1);
            }
            if ( this->_ppacc->_ssss_correlations ) {
              pcps::xscal(_nsites*_nsites*_nsites*_nsites, pcps::zero<SCALAR>(), &_ppacc->_process_ssss.at(0),      1);
            }
            if ( this->_ppacc->_spin_triple_products ) {
              pcps::xscal(                             20, pcps::zero<SCALAR>(), &_ppacc->_process_stp.at(0),       1);
            }
            if ( this->_ppacc->_spin_triangle_3sz ) {
              pcps::xscal(                             20, pcps::zero<SCALAR>(), &_ppacc->_process_3sz.at(0),       1);
            }

            if ( this->_ppacc->_spin_pair_correlations ) {
              pcps::reduce(comm, _thread_sz,        &_ppacc->_process_sz.at(0),                                _nsites, MPI::SUM, 0);
              pcps::reduce(comm, _thread_sz_dot_sz, &_ppacc->_process_sz_dot_sz.at(0),                 _nsites*_nsites, MPI::SUM, 0);
              pcps::reduce(comm, _thread_s_dot_s,   &_ppacc->_process_s_dot_s.at(0),                   _nsites*_nsites, MPI::SUM, 0);
            }
            if ( this->_ppacc->_ssss_correlations ) {
              pcps::reduce(comm, _thread_ssss,      &_ppacc->_process_ssss.at(0),      _nsites*_nsites*_nsites*_nsites, MPI::SUM, 0);
            }
            if ( this->_ppacc->_spin_triple_products ) {
              pcps::reduce(comm, _thread_stp,       &_ppacc->_process_stp.at(0),                                    20, MPI::SUM, 0);
            }
            if ( this->_ppacc->_spin_triangle_3sz ) {
              pcps::reduce(comm, _thread_3sz,       &_ppacc->_process_3sz.at(0),                                    20, MPI::SUM, 0);
            }

            // compute diagonal terms of ssss correlations
            if ( this->_ppacc->_ssss_correlations ) {
              const double s = 0.5 * double(_configs_per_site-1);
              for (int i = 0; i < _nsites; i++)
              for (int k = 0; k < _nsites; k++)
              for (int l = 0; l < _nsites; l++) {
                _ppacc->_process_ssss[_nsites*(_nsites*(_nsites*i+i)+k)+l] = s * ( s + 1.0 ) * _ppacc->_process_s_dot_s[k*_nsites+l];
                _ppacc->_process_ssss[_nsites*(_nsites*(_nsites*k+l)+i)+i] = s * ( s + 1.0 ) * _ppacc->_process_s_dot_s[k*_nsites+l];
              }
            }

          }

          // wait for thread 0
          pcps::pthread_barrier_wait(this->_pacc->barrier());

        }

      }

      // function to compute the local configuration's contribution to the spin triangles' triple products
      void accumulate_spin_triple_products(const long int count, const int * const config,
                                           pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf) {

        // check that the particles are spins
        if (typeid(PART) != typeid(pcps::Part_Spin))
          throw pcps::Exception("accumulate_spin_triple_products only implemented for spins");

        // initialize the work configuration
        for (int i = 0; i < _nsites; i++)
          _workspace[i] = config[i];

        // function object for getting the triple product contributions
        struct Cont_ijk {

          static SCALAR get(const int configs_per_site, pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf,
                            const pcps::action_value_object<SCALAR, pcps::OpRL<SCALAR>, PART> & avoRL,
                            const int * const config, int * const work, const int i, const int j, const int k) {

            // create an array to hold operators' sites' occupations
            int op_occs[2];

            // create an array to hold affected site indices
            int affected[4];

            // get the spin
            const double s = 0.5 * double(configs_per_site-1);

            // initialize the contribution
            SCALAR cont = pcps::zero<SCALAR>();

            // record which sites are changing
            affected[0] = k;
            affected[1] = j;

            // get the contribition from Sk+Sj-
            work[k] -= 1;
            work[j] += 1;
            if ( work[k] >= 0 && work[j] < configs_per_site ) {
              op_occs[0] = work[k];
              op_occs[1] = work[j];
              cont += cvf.get_coeff_ratio(2, affected, work, config) * avoRL.get_value(op_occs, configs_per_site);
            }

            // get the contribition from Sj+Sk-
            work[k] += 2;
            work[j] -= 2;
            if ( work[j] >= 0 && work[k] < configs_per_site ) {
              op_occs[0] = work[j];
              op_occs[1] = work[k];
              cont -= cvf.get_coeff_ratio(2, affected, work, config) * avoRL.get_value(op_occs, configs_per_site);
            }

            // restore work configuration
            work[k] -= 1;
            work[j] += 1;

            // get the Sz contribution
            cont *= -0.5 * pcps::imaginary_unity<SCALAR>() * ( double(config[i]) - s );

            // return the contribution
            return cont;

          }

        };

        // compute the contribution to each triangles' spin triple product
        for (int t = 0; t < this->_ppacc->_triangles.size(); t++) {

          // get the triangle's site indices
          int i, j, k;
          this->_ppacc->abc_order_ico_triangle_sites(this->_ppacc->_triangles.at(t).at(0),
                                                     this->_ppacc->_triangles.at(t).at(1),
                                                     this->_ppacc->_triangles.at(t).at(2), i, j, k);

          // get the three contributions
          _thread_stp[t] += double(count) * Cont_ijk::get(_configs_per_site, cvf, _ppacc->_avoRL, config, _workspace, i, j, k);
          _thread_stp[t] += double(count) * Cont_ijk::get(_configs_per_site, cvf, _ppacc->_avoRL, config, _workspace, j, k, i);
          _thread_stp[t] += double(count) * Cont_ijk::get(_configs_per_site, cvf, _ppacc->_avoRL, config, _workspace, k, i, j);

        }

      }

      // function to compute the local configuration's contribution to the spin triangles' sz*sz*sz
      void accumulate_spin_triangle_3sz(const long int count, const int * const config,
                                        pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf) {

        // check that the particles are spins
        if (typeid(PART) != typeid(pcps::Part_Spin))
          throw pcps::Exception("accumulate_spin_triangle_3sz only implemented for spins");

        // get the spin
        const double s = 0.5 * double(_configs_per_site-1);

        // compute the contribution to each triangles' sz*sz*sz product
        for (int t = 0; t < this->_ppacc->_triangles.size(); t++) {

          // get the triangle's site indices
          const int i = this->_ppacc->_triangles.at(t).at(0);
          const int j = this->_ppacc->_triangles.at(t).at(1);
          const int k = this->_ppacc->_triangles.at(t).at(2);

          // get the contribution
          _thread_3sz[t] += double(count) * pcps::unity<SCALAR>() * ( double(config[i]) - s )
                                                                  * ( double(config[j]) - s )
                                                                  * ( double(config[k]) - s );

        }

      }

  };

  //------------------------------------------------------------------------------------------
  // pcps::PropertiesProcessAccum::create_thread_accum -- Function to create a thread accumulator
  //                                                      from a process accumulator.
  //------------------------------------------------------------------------------------------

  template <class SCALAR, class PART, class REF>
  pcps::ThreadAccum<SCALAR, PART, REF> * pcps::PropertiesProcessAccum<SCALAR, PART, REF>::create_thread_accum(const int tid, pcps::Arena & arena) {
    return new(arena) pcps::PropertiesThreadAccum<SCALAR, PART, REF>(this, tid, arena);
  }

} // end namespace pcps

#endif
